MacHack 1998/MacHack 1998.toast/The Hacks!/PCA Icon Arranger ƒ/PCoord3.p

Text File | 1998-06-20 | 12.8 KB | 538 lines | [TEXT/PJMM]

unit PCoord3;
{ Modification History }
{ March 1998 Think Pascal version by Philippe Casgrain }
{ June 1998 made in a unit for MacHack }

interface

    uses
        Printing, MacInit, Statis, MathDeclarations;

    function GetCoords (resemblanceMatrix: TMatrixPtr;
                                    estUneMatriceDistance: boolean;
                                    nombreDimensions: Integer): TMatrixPtr;

implementation

{    Printing, MacInit, Statis;}
    const
        nbcar = 10;
        MenuCoVAR = 6002;
    type
        VARiint = 1..3;
        VARiable = record
                case VARiint of
                    1: (
                            int: INTEGER
                    );
                    2: (
                            re: ExtENDed
                    );
                    3: (
                            car: alpha
                    );
            end;
    var
        i, j, k, p, nbdesc: INTEGER;
        ptrmin, ptrmax, compte, nbvalides, nbespaces, nobj, ASauter: INTEGER;
        x1, plusgrand, pluspetit, y, y1, rayoncarre, valmax, MaxPt: ExtENDed;
        distance: BOOLEAN;
        VARi: VARiable;
        entree: FileTYPE;
        GLin, GMult, x, b, c, cs, sn, pp, Diagonale, Ptri: Ptr;
        car: char;
        FichPos, FichVP: FileTYPE;
        ThePtr: PtrTYPE;
        grow, MatDim, Space: Size;
        ShowGeneral: WindowPtr;
        Menu2Hdl: MenuHandle;

        resultMatrix: TMatrixPtr;


    procedure symmetric (n: INTEGER;
                                    g, x, a: Ptr);
        label
            11, 77;
        var
            i, j, k, m, m1: LongInt;
            Ancienm: INTEGER;
            t, norm, eps, sine, cosine, lambda, mu, a0, a1, b0, beta, x0, x1, w1, w2, Temp1, Temp2: ExtENDed;

        function max (a, b: ExtENDed): ExtENDed;
        begin
            if a > b then
                max := a
            else
                max := b;
        end;

        procedure householder (n: INTEGER;
                                        var g, x, a, b: Ptr;
                                        var norm: ExtENDed);
            var
                i, j, k: LongInt;
                t, sigma, alpha, beta, gamma, absb, provi, Temp1, Temp2: ExtENDed;
        begin
            norm := 0;
            absb := 0;
            for k := 1 to n - 2 do begin
                with AdVec(a, k).PtrRee^ do
                    Reel := AdLin(G, k, k).PtrRee^.Reel;
                sigma := 0;
                for i := k + 1 to n do
                    with AdLin(g, i, k).PtrRee^ do
                        sigma := sigma + sqr(Reel);
                with AdVec(a, k).PtrRee^ do
                    t := absb + abs(Reel);
                absb := sqrt(sigma);
                norm := max(norm, t + absb);
                with AdLin(g, k + 1, k).PtrRee^ do
                    alpha := Reel;
                if alpha < 0 then
                    beta := absb
                else
                    beta := -absb;
                with AdVec(b, k).PtrRee^ do
                    Reel := beta;
                if sigma <> 0 then begin
                    gamma := 1 / (sigma - alpha * beta);
                    with AdLin(G, k + 1, k).PtrRee^ do
                        Reel := alpha - beta;
                    for i := k + 1 to n do begin
                        provi := 0;
                        for j := k + 1 to i do begin
                            with AdLin(GLin, i, j).PtrRee^ do
                                w1 := Reel;
                            with AdLin(G, j, k).PtrRee^ do
                                provi := provi + w1 * Reel;
                        end;
                        for j := i + 1 to n do begin
                            with AdLin(G, j, i).PtrRee^ do
                                w1 := Reel;
                            with AdLin(G, j, k).PtrRee^ do
                                provi := provi + w1 * Reel;
                        end;
                        with AdVec(pp, i).PtrRee^ do
                            Reel := provi * gamma;
                    end;
                    t := 0;
                    for i := k + 1 to n do begin
                        with AdLin(G, i, k).PtrRee^ do
                            temp1 := Reel;
                        with AdVec(pp, i).PtrRee^ do
                            t := t + temp1 * Reel;
                    end;
                    t := t * 0.5 * gamma;
                    for i := k + 1 to n do begin
                        with AdLin(G, i, k).PtrRee^ do
                            temp1 := Reel;
                        with AdVec(pp, i).PtrRee^ do
                            Reel := Reel - t * temp1;
                    end;
                    for i := k + 1 to n do
                        for j := k + 1 to i do begin
                            with AdLin(G, i, k).PtrRee^ do
                                w1 := Reel;
                            with AdLin(G, j, k).PtrRee^ do
                                w2 := Reel;
                            with AdVec(pp, j).PtrRee^ do
                                Temp1 := Reel;
                            with AdVec(pp, i).PtrRee^ do
                                Temp2 := Reel;
                            with AdLin(G, i, j).PtrRee^ do
                                Reel := Reel - w1 * Temp1 - Temp2 * w2;
                        end;
                    for i := 2 to n do begin
                        provi := 0;
                        for j := k + 1 to n do begin
                            with AdLin(G, j, k).PtrRee^ do
                                W1 := Reel;
                            with AdMat(x, i, j).PtrRee^ do
                                provi := provi + w1 * Reel;
                        end;
                        with AdVec(pp, i).PtrRee^ do
                            Reel := gamma * provi;
                    end;
                    for i := 2 to n do
                        for j := k + 1 to n do begin
                            with AdLin(G, j, k).PtrRee^ do
                                w1 := Reel;
                            with AdVec(pp, i).PtrRee^ do
                                Temp1 := Reel;
                            with AdMat(x, i, j).PtrRee^ do
                                Reel := Reel - Temp1 * w1;
                        end;
                end;
            end;
            with AdLin(G, n - 1, n - 1).PtrRee^ do
                Temp1 := Reel;
            with AdLin(G, n, n).PtrRee^ do
                Temp2 := Reel;
            with AdLin(G, n, n - 1).PtrRee^ do
                w1 := Reel;
            with AdVec(a, n - 1).PtrRee^ do
                Reel := Temp1;
            with AdVec(a, n).PtrRee^ do
                Reel := Temp2;
            with AdVec(b, n - 1).PtrRee^ do
                Reel := w1;
            t := abs(w1);
            norm := max(norm, absb + abs(Temp1) + t);
            norm := max(norm, t + abs(Temp2));
        end;        (* fin de householder *)

    begin        (* debut de symmetric *)
        for i := 1 to n do begin
            with AdMat(x, i, i).PtrRee^ do
                Reel := 1;
            for j := i + 1 to n do begin
                with AdMat(x, i, j).PtrRee^ do
                    Reel := 0;
                with AdMat(x, j, i).PtrRee^ do
                    Reel := 0;
            end;
        end;
        with AdVec(b, n).PtrRee^ do
            Reel := 0;
        householder(n, g, x, a, b, norm);
        eps := norm * 1.3e-7;
        with AdVec(b, 0).PtrRee^ do
            Reel := 0;
        mu := 0;
        m := n;
        AncienM := 0;
11:
        if m = 0 then
            goto 77
        else begin
            i := m - 1;
            k := m - 1;
            m1 := m - 1;
        end;
        if AncienM <> m then begin
            AncienM := M;
        end;
        with AdVec(b, k).PtrRee^ do
            w1 := Reel;
        if abs(W1) <= eps then begin
            with AdVec(a, m).PtrRee^ do
                w1 := Reel;
            with AdLin(G, m, m).PtrRee^ do
                Reel := w1;
            m := k;
            goto 11;
        end;
        repeat
            i := i - 1;
            with AdVec(b, i).PtrRee^ do
                w1 := Reel;
        until abs(w1) <= eps;
        k := i + 1;
        with AdVec(b, m1).PtrRee^ do
            b0 := sqr(Reel);
        with AdVec(a, m1).PtrRee^ do
            Temp1 := Reel;
        with AdVec(a, m).PtrRee^ do
            Temp2 := Reel;
        a1 := sqrt(sqr(Temp1 - Temp2) + 4 * b0);
        t := Temp1 * Temp2 - b0;
        a0 := Temp1 + Temp2;
        if a0 >= 0 then
            lambda := a0 + a1
        else
            lambda := a0 - a1;
        t := t / lambda;
        if abs(t - mu) < 0.5 * abs(t) then begin
            mu := t;
            lambda := t;
        end
        else if abs(lambda - mu) < 0.5 * abs(lambda) then
            mu := lambda
        else begin
            mu := t;
            lambda := 0;
        end;
        with AdVec(a, k).PtrRee^ do
            Reel := Reel - lambda;
        with AdVec(b, k).PtrRee^ do
            beta := Reel;
        for j := k to m1 do begin
            with AdVec(a, j).PtrRee^ do
                a0 := Reel;
            with AdVec(a, j + 1).PtrRee^ do
                a1 := Reel - lambda;
            with AdVec(b, j).PtrRee^ do
                b0 := Reel;
            t := sqrt(sqr(a0) + sqr(beta));
            cosine := a0 / t;
            with AdVec(cs, j).PtrRee^ do
                Reel := cosine;
            sine := beta / t;
            with AdVec(sn, j).PtrRee^ do
                Reel := sine;
            with AdVec(a, j).PtrRee^ do
                Reel := cosine * a0 + sine * beta;
            with AdVec(a, j + 1).PtrRee^ do
                Reel := -sine * b0 + cosine * a1;
            with AdVec(b, j).PtrRee^ do
                Reel := cosine * b0 + sine * a1;
            with AdVec(b, j + 1).PtrRee^ do
                beta := Reel;
            with AdVec(b, j + 1).PtrRee^ do
                Reel := cosine * beta;
            with AdVec(c, j).PtrRee^ do
                Reel := sine * beta;
        end;
        with AdVec(b, k - 1).PtrRee^ do
            Reel := 0;
        with AdVec(c, k - 1).PtrRee^ do
            Reel := 0;
        for j := k to m1 do begin
            with AdVec(sn, j).PtrRee^ do
                sine := Reel;
            with AdVec(cs, j).PtrRee^ do
                cosine := Reel;
            with AdVec(a, j).PtrRee^ do
                a0 := Reel;
            with AdVec(b, j).PtrRee^ do
                b0 := Reel;
            with AdVec(c, j - 1).PtrRee^ do
                w1 := Reel;
            with AdVec(b, j - 1).PtrRee^ do
                Reel := Reel * cosine + w1 * sine;
            with AdVec(a, j).PtrRee^ do
                Reel := a0 * cosine + b0 * sine + lambda;
            with AdVec(b, j).PtrRee^ do
                Reel := -a0 * sine + b0 * cosine;
            with AdVec(a, j + 1).PtrRee^ do
                Reel := Reel * cosine;
            for i := 1 to n do begin
                with AdMat(x, i, j + 1).PtrRee^ do
                    x1 := Reel;
                with AdMat(x, i, j).PtrRee^ do begin
                    x0 := Reel;
                    Reel := x0 * cosine + x1 * sine;
                end;
                with AdMat(x, i, j + 1).PtrRee^ do
                    Reel := -x0 * sine + x1 * cosine;
            end;
        end;
        with AdVec(a, m).PtrRee^ do
            Reel := Reel + lambda;
        goto 11;
77:
    end;  (* Fin Symmetric *)


    procedure nouvellescoord;
        var
            i, j, k: INTEGER;

    begin
        if NewMatrix(resultMatrix, nbDesc, nbEspaces, false, false) then
            for j := 1 to nbdesc do begin
                for i := 1 to nbespaces do begin
                    with AdVec(Ptri, i).PtrInt^ do
                        k := Int;
                    with AdMat(x, j, k).PtrRee^ do begin
                        SetElement(resultMatrix, j, i, Reel);
                    end;
                end;
            end;
    end;        (* fin nouvelles coordonnees *)




    procedure shell (a, Ptri: Ptr;
                                    n: INTEGER);
        label
            1;
        var
            i, j, k, m, o, i1, i2: INTEGER;
            w, ww: ExtENDed;
    begin
        for i := 1 to n do
            m := 2 * i - 1;
        m := m div 2;
        while m <> 0 do begin
            k := n - m;
            for j := 1 to k do begin
                i := j;
                while i >= 1 do begin
                    with AdVec(a, i).PtrRee^ do
                        w := Reel;
                    with AdVec(a, i + m).PtrRee^ do
                        if Reel <= w then
                            goto 1;
                    o := i + m;
                    with AdVec(a, o).PtrRee^ do
                        ww := Reel;
                    with AdVec(a, i).PtrRee^ do
                        Reel := ww;
                    with AdVec(a, o).PtrRee^ do
                        Reel := w;
                    with AdVec(Ptri, o).PtrInt^ do
                        i2 := Int;
                    with AdVec(Ptri, i).PtrInt^ do
                        i1 := Int;
                    with AdVec(Ptri, i).PtrInt^ do
                        Int := i2;
                    with AdVec(Ptri, i).PtrInt^ do
                        Int := i1;
                    i := i - m;
                end;
1:
            end;
            m := m div 2;
        end;
    end;    (* fin de shell *)

    procedure Normalise;
        var
            i, j, k: INTEGER;
            w: ExtENDed;
    begin
        for j := 1 to nbdesc do begin
            for i := 1 to nbValides do begin
                with AdVec(Diagonale, i).PtrRee^ do
                    w := Reel;
                with AdVec(Ptri, i).PtrInt^ do
                    k := Int;
                with AdMat(x, j, k).PtrRee^ do
                    Reel := Reel * sqrt(w);
            end;
        end;
    end;

    procedure CopyVP;
        var
            i: INTEGER;
    begin
        x1 := 0;
        y := 0;
        for i := 1 to nbdesc do
            with AdVec(Diagonale, i).PtrRee^ do
                x1 := x1 + Reel;
        with AdVec(Diagonale, NbDesc).PtrRee^ do
            if Reel < 0 then begin
                y := -Reel;
                x1 := x1 + nbdesc * y;
            end;
    end;


    procedure EcritEntete (matRessemblance: TMatrixPtr);
        var
            I, J, k: LongInt;
            Mult: Variable;
            x, Total, w, ww: Extended;

    begin
        Nbdesc := matRessemblance^.n;
        NobjSimil := NbDesc;
        Space := NbDesc;
        Space := (Space * (Space + 1)) div 2;
        MatDim := (NbDesc + 1) * 3;
        if Space < MatDim then
            Space := MatDim;
        Diagonale := Memoire(1, NbDesc, 1, 1, ReelBytes, True);
        if NbDesc mod 2 = 0 then
            GLin := Memoire(1, NbDesc div 2, 1, NbDesc + 1, ReelBytes, True)
        else
            GLin := Memoire(1, NbDesc, 1, (NbDesc + 1) div 2, ReelBytes, True);
        Nobj := matRessemblance^.n; { pas important }
        NbDescSimil := Nobj;

        for i := 1 to nbdesc do
            with AdVec(Diagonale, i).PtrRee^ do
                Reel := 0;
        total := 0;
        k := 0;
        (* remplissage de la matrice *)
        for i := 1 to nbdesc do begin
            for j := i + 1 to nbdesc do begin
                Mult.re := GetElement(matRessemblance, i, j);
                if distance then
                    x := (-sqr(Mult.re) / 2)
                else
                    x := (-sqr(1 - Mult.re) / 2);
                with AdLin(GLin, j, i).PtrRee^ do
                    Reel := x;
                with AdVec(Diagonale, i).PtrRee^ do
                    Reel := Reel + x;
                with AdVec(Diagonale, j).PtrRee^ do
                    Reel := Reel + x;
                total := total + x + x;
            end;
            with AdLin(GLin, i, i).PtrRee^ do
                Reel := 0;
        end;
        for i := 1 to nbdesc do
            with AdVec(Diagonale, i).PtrRee^ do
                Reel := Reel / nbdesc;
        total := total / nbdesc / nbdesc;
        for i := 1 to nbdesc do begin
            for j := 1 to i do begin
                with AdVec(Diagonale, i).PtrRee^ do
                    w := Reel;
                with AdVec(Diagonale, j).PtrRee^ do
                    ww := Reel;
                with AdLin(GLin, i, j).PtrRee^ do
                    Reel := Reel - w - ww + total;
            end;
        end;
    end;    (* ENTETE *)

    function GetCoords (resemblanceMatrix: TMatrixPtr;
                                    estUneMatriceDistance: boolean;
                                    nombreDimensions: Integer): TMatrixPtr;
        var
            fatalError: Boolean;
            i: Integer;
    begin    (* prog princ *)
        fatalError := false;
        resultMatrix := nil;

        FichierSimil := TRUE;
        NbFiles := 0;
        Distance := estUneMatriceDistance;
{LisFichSimil(4, entree, TRUE); voila qui ouvre une ref au fichier simil dans entree}
        EcritEntete(resemblanceMatrix);
        x := Memoire(1, NbDesc, 1, NbDesc, ReelBytes, TRUE);
        b := Memoire(0, NbDesc + 1, 1, 1, ReelBytes, TRUE);
        c := Memoire(0, NbDesc + 1, 1, 1, ReelBytes, TRUE);
        cs := Memoire(0, NbDesc + 1, 1, 1, ReelBytes, TRUE);
        sn := Memoire(0, NbDesc + 1, 1, 1, ReelBytes, TRUE);
        pp := Memoire(0, NbDesc + 1, 1, 1, ReelBytes, TRUE);
        symmetric(nbdesc, glin, x, diagonale);
        DisposeMemoire(b);
        DisposeMemoire(c);
        DisposeMemoire(cs);
        DisposeMemoire(sn);
        DisposeMemoire(pp);
        GMult := Memoire(1, NbDesc + 2, 1, 3, ReelBytes, TRUE);
        Ptri := Memoire(1, NbDesc, 1, 1, IntBytes, TRUE);
        for i := 1 to nbdesc do
            with AdVec(Ptri, i).PtrInt^ do
                Int := i;
        shell(diagonale, Ptri, nbdesc);
        nbvalides := 0;
        while (AdVec(Diagonale, nbvalides + 1).PtrRee^.Reel > epsilon) and (nbvalides < nbdesc) do
            nbvalides := nbvalides + 1;
        Normalise;
        CopyVP;
        if nbvalides <= 1 then
            fatalError := true;
        if not fatalError then begin
            nbespaces := nombreDimensions;
            if NbEspaces > 0 then begin
                NouvellesCoord;
            end;
        end; { if not fataError }

        GetCoords := resultMatrix;
    end; { GetCoordinates }

end. { unit }